function fv = calib_err2(x)
% input: parameters to be calibrated including pro, z0, c
% output: errors of moments to be matched in line 67 fv = sqrt(sum(fvv(1:nn).^2)); 

pro = x(1); % total factor productivity
z0 = x(2); % productivity of new entrants
c = x(3); % entry costs
% nn = 4 means matching bank size, nn = 3 not considering bank size

%rho = x(4);
global beta mu_z sig_z lambda I m_bar eta z_min z_max rho nn

R = rate(pro,z0,c);

z = linspace(z_min,z_max,I)';
dz = z(2) - z(1);
dz2 = dz^2;
i_entry = max(round((z0 - z_min)/dz),1);
z0 = z(i_entry);
psi = zeros(I,1);
psi(i_entry) = 1/dz;

qstar = z + beta*R - 1;
piz = exp(qstar)/beta;
    
    %% HJB
mu = mu_z*ones(I,1);
sig = sig_z*ones(I,1);
sig2 = sig.^2;
    
X = -min(mu,0)/dz + sig2/(2*dz2);
Y = -max(mu,0)/dz + min(mu,0)/dz - sig2/dz2;
Z = max(mu,0)/dz + sig2/(2*dz2);
A = spdiags(Y,0,I,I) + spdiags(X(2:I),-1,I,I) + spdiags([0;Z(1:I-1)],1,I,I);
    
A(1,1) = Y(1) + X(1); % Reflecting barrier v'(z_min) = 0
A(I,I) = Y(I) + Z(I); % Reflecting barrier v'(z_max) = 0
    
B = (rho+lambda)*speye(I) - A;
v = B\piz; % value functions derived from HJB
    
    %% KFE
    % entry
m = m_bar*exp(eta*(sum(v.*psi.*dz)-c));
%m = m_bar; % reach steady state

Aa = A';
% Aa(1,1) = 0; % Reflecting barrier f(z_min) = 0
% Aa(I,I) = 0; % Reflecting barrier f(z_max) = 0
Bb = lambda*speye(I) - Aa;
f = Bb\(m.*psi); % stationary distribution of z derived from KFE

%%
entry_size = exp(z0+beta*R-1);
% bank_size = sum(exp(qstar).*f.*dz)/sum(f.*dz);
% market_to_book = sum(v./exp(qstar).*f.*dz)/sum(f.*dz);
market_to_book = sum(v.*f.*dz)/sum(.1*exp(qstar).*f.*dz); % weighted average, 10 times leverage

% fv(1) = bank_size - 11;
% fv(1) = (market_to_book/1.1-1);
fvv = zeros(6,1);
fvv(1) = R/0.05-1; %lending rate
fvv(2) = entry_size/0.25-1; % size of entrants
fvv(3) = (market_to_book/1.1-1);
fvv(4) = sum(exp(qstar).*f.*dz)/sum(f.*dz)/19.067-1; % average bank size
%fvv(5) = log(250)+1-beta*R - z_max; 
%fvv(6) = log(0.1)+1-beta*R - z_min;

fv = sqrt(sum(fvv(1:nn).^2));

disp(['Total factor productivity: ', num2str(pro)])
disp(['Entry productivity: ', num2str(z0)])

disp(['size of entrants: ', num2str(exp(z0+beta*R-1))])
disp(['average bank size: ', num2str(sum(exp(qstar).*f.*dz)/sum(f.*dz))])
disp(['average market to book: ', num2str( market_to_book)])
disp(['mass: ',num2str(sum(f))])
disp(['Upper bound of bank size: ', num2str(max(exp(qstar)))])
disp(['Lower bound of bank size: ' ,num2str(min(exp(qstar)))])
disp(['Entry cost: ',num2str(sum(v.*psi.*dz))])
disp(['Lending rates: ', num2str(R)])